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Abstract 

We here give an improved formalism for calculating the evolution of density fluctuations 
and temperature perturbations in flat universes. Our equations are general enough to treat 
the perturbations in collisionless relics like massive neutrinos. We find this formulation to be 
simpler to use than gauge dependent and other gauge-invariant formalisms. We show how 
to calculate temperature fluctuations (including multipole moments) and transfer functions, 
including the case of collisionless relics like massive neutrinos. We call this formalism “quasi- 
Newtonian” because the equations for the potential and cold matter fluctuation evolution have 
the same form as the Newtonian gravitational equations in an expanding space. The density 
fluctuation variable also has the same form inside and outside of the horizon which allows the 
initial conditions to be specihed in a simple intuitive way. Our sample calculations demonstrate 
how to use these equations in cosmological models which have hot, cold, and mixed dark matter 
and adiabatic (isentropic) or isocurvature modes. We also give an approximation may be used 
to get transfer functions quickly. 

Subject Headings: cosmology: theory - dark matter - cosmic microwave background 



1 Introduction 


Finding the correct model for the formation of structure in the universe is currently a subject of 
intense study. It is generally believed that cosmic structures are the result of the gravitational 
amplification of initially small density fluctuations in the early universe. There are candidate 
mechanisms for producing the primordial fluctuations, e.g., inflation and topological defects, which 
must be tested against observations of the current universe. To test these primordial mechanisms, 
one has to evolve the perturbations up to the present day. This task is complicated by the fact 
that the amount and composition of dark matter is not known, so potentially a large number of 
possibilities may exist which need to be checked against observations. 

The two main pieces of information required for testing theories are the predicted power 
spectra of the temperature and matter fluctuations. Various prescriptions for how to calculate 
these exist in the literature (Kodama & Sasaki, 1984; Sasaki & Gouda, 1986; Bond & Efstathiou, 
1987; Mukhanov, Feldman, & Brandenberger, 1987; Durrer &: Straumann, 1988; Stompor, 1994). 
The problem with these treatments is that either they do not span the current range of models or 
they that use variables which are unsatisfying in that the physical nature of the initial conditions 
in is not transparent. 

Many of the treatments do not include equations for how to deal with the free streaming of 
massive collisionless relics such as few eV mass neutrinos, an element which is indicated in critical 
density universe models based on inflation (Schaefer & Shah, 1994; Pogosyan & Starobinsky, 1994; 
Fiddle & Lyth, 1994). Some gauge-invariant treatments which include a treatment of collisionless 
relics (Durrer & Straumann, 1988; Stompor, 1994) use a generalized density perturbation vari¬ 
able (there are several reasonable candidates) which has a complicated behavior outside of the 
horizon and thus makes it inconvenient for specifying initial conditions. Here we present a unified 
treatment of the evolution of matter and temperature perturbations, which include free streaming 
massive neutrinos, and use a density perturbation variable which is most natural to use from the 
point of view of matter perturbations. We also have attempted to put this into a form which can 
be easily used by someone to get results without the need to slog through the subtleties of the 
general relativistic underpinnings of the formalism. 

For the sake of brevity we have only treated flat universes here. We expect to follow this paper 
with another which covers universes with curvature. 

We begin in section 2 by constructing a gauge-invariant definition of the distribution function, 
in a slightly different form than that of Durrer &: Straumann, (1988). We then recast the generally 
covariant Boltzmann equation for this distribution function into a more useful form, and this is 
the basis for the rest of the paper. We also derive fluid equations for the first three moments of 
the distribution function. In section 3 we discuss how to treat relativistic collisionless particles, 
e.g. massless neutrinos and decoupled photons, using only fluid variable equations. Section 4 
explains how to derive initial conditions. In section 5 we discuss how to treat the baryon-photon 
fluids which are coupled early on, and then show how to calculate the present day temperature 
anisotropies. Section 6 gives examples of numerical implementation of these equations, and some 
details about the numerical computations. We end with a summary of our main results. Lastly, 
there is an appendix which gives a quick way to calculate transfer functions for some models. 
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2 General Gauge-Invariant Evolution Equations for Gollisionless 
Relics 

Collisionless massive relic particles cannot be adequately described in a fluid formulation of the 
gravitation equations, (unless the particles are “cold” on the scales of interest, i.e. their thermal 
energy is much less than the gravitational potential of a density perturbation. One must go back 
to the Boltzmann equation. Since we are working in a regime in which general relativistic effects 
are important, the best place to start is with the generally covariant Boltzmann equation. This 
equation has been derived and discussed in many places (see e.g., Lindquist, 1966; Ehlers, 1969; 
Stewart, 1971; Kodama & Sasaki, 1984; Durrer & Straumann, 1988), and we need not repeat the 
details here. [An interesting approach to solving this equation in power series form has also been 
explored (Rebhan & Schwartz, 1994).] We will only sketch the pertinent parts of the derivation 
of the perturbed Boltzmann equation. Small perturbations of the distribution function of the 
particles are not in general gauge invariant. We will remove the gauge dependent part from 
the distribution function. Next we expand the gauge-invariant distribution function in terms of 
angular moments, the first three of which are associated with gauge invariant fluid variables: 
the density, velocity, entropy, and anisotropic pressure perturbations. We will then discuss the 
correspondence of the initial values of these moments to the initial fluid variable values. We end 
this section with some remarks on the physical interpretation of the equations. 

2.1 The Generally Covariant Boltzmann Eqnation 

The non~relativistic distribution function specifies the number of particles at time t with velocity 
between v and v^dv and position between x and x^dx. Since we intend to describe the behavior 
of relativistic particles in a perturbed spacetime, we need to use a generally covariant formulation, 
so we must take care that our definition of the distribution function is also generally covariant. 
For the setup of the Boltzmann equation, we follow the treatment of Kodama & Sasaki, (1984). 
We will also borrow from Durrer & Straumann, (1988). 

In the general relativistic version of kinetic theory, we use the invariant volume elements to 
define the distribution function. Here we will need the space-like covariant version of the volume 
element a a'. 

a a = --^UaU^^e^^Xadx’^dx^dx'' (1) 

where is an arbitrary time like vector field and is the totally antisymmetric tensor with 

eoi 23 = \/(—det^^i,). is the metric. Greek indices run from 0 to 3, while latin indices i through 
m will run from 1 to 3. The other latin indices will be used to distinguish the different component 
fluids. The invariant momentum space volume element for a particle with mass m is TTq 

TTg = + m^)e^yXadq^dq'' dq^dq"^. (2) 

Here q^ is the contravariant momentum, the 5 function enforces the relation among energy mass 
and momentum (keeps the particle momenta “on the mass shell”), and the Heavyside {9) function 
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restricts us to the regime of positive energies for the particles. The number of particles crossing 
a unit hypersurface element cJq, within a momentum space volume tt^ is given by 


dn = /(x", q^)q^crf,TTg. (3) 

Conservation of q^cr^TTq along a particle trajectory is guaranteed by the Liouville theorem. Any 
change in dn due to interparticle collisions can be represented by 

C{f)q^crf,TTg. (4) 

The relativistic Boltzmann equation becomes 


^(/) = C{f). 


(5) 


where C{f) is the linear operator which is the total derivative of the distribution function with 
respect to an affine parameter A which traces geodesics. In non-relativistic theory A would be 
identified with the time variable. In our relativistic version, C is 


d ,, d dq>^ d 

— = q^ -h — - 

dX dx>^ dX dq>^ 


( 6 ) 


where the momentum is defined, as usual, by 


qf^ = 


dx^ 

~dX' 


and the particles follow geodesics which satisfy 


dq^^ 

lix 




(7) 

( 8 ) 


The are the metric connection coefficients. Since we are initially interested in collisionless 
particles we will set C{f) = 0. In section 5 we will revisit the collision term for baryon-photon 
scattering. 

Eventually, we will need to relate this microscopic distribution function to macroscopic quan¬ 
tities like the density and pressure. This is done by taking momentum averages. The energy- 
momentum tensor can be obtained from the formula 


= I 


(9) 


In order to evaluate these integrals and solve the Boltzmann equation, one usually introduces a 
tetrad (or vierbein) frame which is convenient for the purpose, usually the one which looks like 
flat space so that dot products are vectors contracted with the Minkowski metric. Kodama & 
Sasaki (1984) introduce a set of tetrads on which the momentum is physical and redshifts with the 
universal expansion. Bond & Szalay (1983) and Durrer &: Straumann (1988) use a tetrad frame 
on which the momentum are comoving in flat expanding space, which has the advantage that the 
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spatial momentum components are constant in time. Here we will follow the latter convention, 
and refer interested readers to the original works for more details about the Boltzmann equation. 

One can define in a covariant way the one-particle distribution function / as a scalar function 
of (conformal) time r, the comoving three-space coordinates x*, and the three-momentum p* (for 
particles on the mass shell). We want to consider perturbations on a homogeneous and isotropic 
background universe, so we split the distribution function into background -|- perturbation. The 
spatial dependence of the perturbations can be expanded in terms of harmonic functions Y (/c*, x*), 
which for flat space (e.g., critical density universes) is like Fourier transforming the equations, as 
Y = We then split the distribution function into 

f = f + i5f)Y (10) 

where the perturbation 6f depends on the Fourier wavevector /c*, the time r and the three 
momentum. We also perturb the metric as 

~ ® T h^i/^ (If) 


where a is the time dependent “scale factor”, and is the Minkowski metric. The scale factor 
satisfies the equation 

a\^ SttG 2 Aa^ 

Qi J 3 3 

where G is Newton’s gravitation constant, p is the total density of mass in the universe, and A 
is the cosmological constant. In the above equation and throughout the paper, the dot (•) over 
a variable denotes a derivative with respect to conformal time r. The conformal time is more 
convenient for doing these calculations than the real time t. The conformal time is, in fact, the 
co-moving horizon size. This property is useful for relating the scale of different effects to the 
scale of structures in the present universe. 

The metric perturbation can also be expanded in terms of harmonic functions: 



^00 

hoj 

hij 


-2A{k,T)Y 

i^B{k,T)Y 

K 

-2HL{k,T)Y5ij 


2HT{k,T) 



kjki 


A:2 


y. 


(13) 


The energy momentum tensor is also perturbed. We can write the perturbations in terms 
of the usual (gauge dependent) fluid perturbations ((5 for the density perturbation, Vf for the fluid 
velocity perturbation, and or tt^ for the isotropic or anisotropic pressure perturbations) and 
T: 


^0 

-^0 


n 


rpl 


-p[l + S{k,T)Y] 
k ■ 

-i^{p + p)vf{k,T)Y 


-P 


t)Y 5j + TTrik, r) 



kjkj \ _ 

k^ ) 


(14) 
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Following Durrer &: Straumann (1988), we will write the Boltzmann equation (in the tetrad 
frame) in terms of the comoving momentum v = p/T, where p is the magnitude of the three vector 
momentum and T is the temperature parameter of the particles. The Boltzmann equation with a 
more physical definition of the momentum which redshifts with time, such as found in (Kodama 
&: Sasaki, 1984; Schaefer, 1991) while being more physically accurate is more complicated to solve 
numerically, so we avoid it. The particle energy variable q (following Durrer & Straumann, 1988) 
associated with the comoving momentum is 


q = 



(15) 


Of course we are using units for which c = = 1. The physical momentum p and the temperature 

T both decrease with time as 1/a, so that v is time independent, but the energy variable is not. 
In this notation the density and pressure are given by 


P 


= T^ J dQdvv'^qf, 


(16) 


and ^ 

P = [ dQdv—f. 

3 7 q 


(17) 


Since on average the universe is homogeneous and isotropic, the background distribution 
function must also be homogeneous and isotropic and is fixed when (and if) the particles were in 
thermal equilibrium in the very early universe: 


f 


1 1 

{2Trh)^ lie’' 


(18) 


where the + is for fermions and the — is for bosons, and h is Planck’s constant divided by 27r. 


2.2 The Gauge—Invariant Perturbed Boltzmann Equation 

The Boltzmann equation for the perturbation to / in eq. dD is then 

((5/)' + ikp-6f = 


q 


df 


ikpqA + v{Hl + -Ht) 

O 


p‘^v{Ht - kB) 


(19) 


The perturbation variables (d/. A, B, and Hl) are gauge dependent, and so are not the most 
physically relevant variables to work with. Durrer & Straumann (1988) showed how to construct 
a gauge invariant version of the distribution function perturbation. The full procedure for doing 
this is given there. Operationally, all we need to know is that we must add something to 5f so 
that the appropriate integrals of 6f yield the gauge-invariant density and velocity perturbations. 
The construction in Durrer & Straumann (1988) was not unique, as one could subtract off the 
gauge dependence in other ways. In fact they chose the combination such that when one converts 
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the distribution function perturbation to its corresponding fluid variable density perturbation by 
integrating over momentum, one would get the density perturbation variable called by Kodama 
& Sasaki, (1984), eg by Bardeen, (1980), and A by Mukhanov, Feldman, &: Brandenberger, (1991). 

A.g =5 + 3(1+ w) + 

where w = p/p. Instead we would like to make another choice for a gauge-invariant distribution 
function in which the corresponding gauge invariant density perturbation is the variable called A 
by Kodama & Sasaki (or Aca for a particular component-note that Kodama & Sasaki also define 
a component perturbation variable called A^, which is not the one we are interested in here.) 



A = 5 + 3(1 + w)-{vf — B) 


( 21 ) 


This variable is in Bardeen’s notation. It is the density perturbation on the spacelike hyper¬ 
surfaces which are the local matter rest frame, and hence is the most natural variable to use from 
the point of the matter fluctuations. Using this variable, the fluid equations most resemble their 
Newtonian counterparts for a fluid under the action of a gravitational perturbation in expanding 
coordinates. Also, this is the density perturbation which Hu & Sugiyama (1995) used in order 
to untangle the various effects which contribute to the cosmic microwave anisotropy. We will 
comment on this aspect in more detail later. We now set about rewriting the Liouville equation 
in a gauge-invariant way, and show explicitly how to use it to calculate the growth of density 
perturbations in critical density universes. 

The gauge invariant generalization T of the distribution function perturbation 5f is given here 


by 




df lav 
dv a k 


{vf - B) + ip^{HT - kB) 

K 


( 22 ) 


Plugging this form into the Boltzmann equation, we get (after some algebraic manipulation) 


B + 


q dv 


df 


V a 


dv 


'3q 

' d' 

2k 

.a. 


1 + w a 
2 


- ciA + rcT-rcH 


a V 

(A + 2tcn) H- V 

a q 


= 0 


(23) 


where all the terms are now manifestly gauge invariant. In the above V is now the gauge- 
invariant fluid velocity perturbation, T is the perturbation in the entropy of the system, and H 
is the amplitude of the anisotropic pressure perturbation. The “sound velocity” is defined 
as Cg = p/p. Note that there is really only one obvious choice for the gauge-invariant velocity 
perturbation V = Vf — Hx/k (Bardeen, 1980; Kodama & Sasaki, 1984; Lyth & Bruni, 1994). The 
other terms T = ttl — 0^6/w and = are already gauge-invariant. 


2.3 Angular Moments of the Distribution Function 

The equation for T can be solved directly, and then one has to do a double integral over p and 
v to get quantities like the density perturbation and velocity perturbation as in (Durrer, 1989). 
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However, one can instead expand the angular dependence of the distribution function perturbation 
in orthogonal polynomials, namely Legendre polynomials (Valdarnini & Bonometto, 1985; 

Schaefer, 1991) as 

-| OO 

= — ^(2£ +l)i^P£(^)crf(A:,r) (24) 

e=o 

and plugging this form into the Boltzmann equation above. Expanding the product + 
l)/iPf(/i) = iPg-i + {i + l)P£+i, resumming the series, and projecting out the coefficients of 
Pi, we get the following set of equations for the cj^ 


di + k- 

q 

— dvr 


£ 


21 + 1 

df 


o'e-i 


dv 


V 


£ + 1 
2£ + l 
2 


o'e+i 


Cl / 9 2 T-r, 

T—-(CgA + iuE - -wH) 

1 -\- w Cl 3 




a. 


(A + 2u;n) + —-E 
q a 


(25) 


where Sij is the Kronecker delta and cj_i = 0. The fluid variables, shown here for component a, 
can be obtained now from integrals over only the magnitude of the momentum v: 


^4 POO 

Aca = — dvv^qao (26) 

pa Jo 


Va = — 


a 


pa Pa J 0 


dvv'^ai 


^4 POO 

WaP-a = — dv -CJ2 

Pa Jo q 


^4 POO 

VJaTa = — dv 

pa Jo 


^ 2 2 


O-Q 


(27) 

(28) 
(29) 


The density perturbation is Aca, the velocity perturbation is 14, Ha is the anisotropic pressure 
perturbation, and Ta is the internal entropy perturbation, which becomes non-zero when density 
and pressure perturbations get out of phase. One can derive from the above equations the values 
for the total perturbations for a multi-component system, finding 


A = PaAca 

(30) 

V = (p + p)"^Xl(^“ 

(31) 

n = p'^^PaHa 

(32) 

r = P y ] I^PaLa + (Ca ~ C^)PaAca 

(33) 
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Using the above relations, the expression in equations (23 & 25) which is common to all the 
component fluid equations can be simply expressed as 


2 2 
lA + wT - -ten = P~^ PaiclAca + WaTa - -rCaHa). 


(34) 


One disadvantage of this formalism is that one must now solve a possibly large number of coupled 
ai equations. However, one can regulate the number of cr^ calculated to get an accurate answer, 
which can be estimated ahead of time. Specifically, one would like to ensure that the logarithmic 
derivative d log cr^/d log r is small. Roughly, this condition is met when fcr(T/m)/t' <C 1 for a 
massive particle (if the particle is massless, just ignore the T/m term) over the entire range of 
integration. For a particle which becomes non-relativistic prior to the epoch of matter/radiation 
equality, one can easily show that the constraint is most stringent when evaluated at equality, 
thus one should satisfy the condition /cr(T/m)/£|eg 1 to ensure an accurate integration. Should 

one be interested in an extremely fast integration in which lower accuracy can be tolerated, one 
can use the abbreviated fluid version of these equations as described in the Appendix. 

One can then integrate the moments [using the procedure as given, for example in Schaefer 
(1991)] to recover the fluid equations given in other work (Bardeen, 1980; Kodama & Sasaki, 
1984). For example, multiplying the equation for (Tq by T^v^q/pa and integrating over v we get 


+ 3-(c„ - Wa)Aca = -3-WaTa 

a a 


— k{l + WajVa + 3--—I- 

a I + w 


czA + rcF — -tcH 


Similarly we obtain from the equation for cji 


Ha + -Va = 2>-cl{Va-V)-^{-) (A + 2u;n) 
a a zk \aJ 


+- 


k 


CqAcjj -|- IVaTa g^aHa 


(1 + Wa) . 

This equation can also be expressed in the following useful form: 

k 


(U-C) + = 


C-nA^a T Wa^a 


k 


cZA + rcF-rcH 

3 


(35) 


(36) 


(37) 


(1 + w) 

For completeness, we also give the equations for the total density and velocity perturbations which 
can be derived using eqs. (30-33, 35, and 


A — 3w-A = —k{l + w)V — 2—wIi^ 
a a 


( 38 ) 





















and 




+ 


k 


(l + w) 


c, A + uiF — -uin 


(39) 


2.4 Initial Values For cr^ 


One can see from the equation (25) for the ai that only the i = 0 and i = 1 modes are driven 
by gravitational coupling to the matter fluctuations in the universe. If we look at the state of the 
free streaming particles when the perturbation is well outside the horizon (kr <C 1), we can see 
that the higher moment equations (pSl) are solved by 


di ~ {krY Vi, 


(40) 


in this regime. However, as we will discuss in section 4 on initial conditions, do ~ krai, so in 
fact the i = 2 mode a 2 ~ do and must be included in the definition of initial conditions. We can 
neglect all I" > 3 moments, however. 

If we start the evolution at r* when all scales of interest are outside the horizon /cr* <C 1, then 
we only need to specify the first three moments. From equation (^) for we see J~{Ti) ought 
to have an initial form like (see also Peebles, 1973) 


df 


—-H + iqiiB + -C'-(3/i^ — 1) 


(41) 


where the factors u/3 and q have been chosen for the simplest interpretation of H, B, and C. 
From the definitions of Aca, 14, and H^ we can see that 


A = 


Ar 


i + Wa 

B = Vain) 


b) 


and 


C = 


5 Wa na(ri) 

3 cl I + Wa 


(42) 

(43) 

(44) 


Using the definition of T in eq. (24) we can now make our identification of the starting values 
for the oi'. 


O-Q 

_ ^df Aca 

3 dvl + Wa 

(45) 

<7l 

4vr df 

(46) 

02 

_ dvr 4/ Wa Ha 

45 dv cl 1 + Wa’ 

(47) 

Oe 

= 0 for I' > 3 

(48) 
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where all variables are to be evaluated at the initial time r^. 

The initial conditions needed for setting up the fluid variables A^a, 14, and Ila will be discussed 
in section ^ We now have all the ingredients necessary for solving the coupled set of equations 
for an arbitrary composition of different matter components for the universe. 

2.5 Why “Quasi-Newtonian”, and the Physical Interpretation of the Equa¬ 
tions 

The main reason for the calling our formalism “quasi-Newtonian” is that the equation for the 
gravitational potential has exactly the Newtonian form. The gauge-invariant gravitational po¬ 
tential <1> and our density perturbation variable are related via Poisson’s equation, which in Fourier 
space is: 

— (J) = dvrGpA (49) 

With other choices for a gauge-invariant density perturbation, this equation only holds on sub¬ 
horizon scales. Since the initial conditions are usually specified when the perturbation scales are 
larger than the horizon, it is convenient for the variables to behave in a simple intuitive manner 
outside the horizon. The equations in these variables then become easily interpretable in terms 
of physics (of the sub-horizon behavior) in this formulation. 

The Newtonian resemblance of this formalism goes deeper than the Poisson relation. The 
equations themselves have the characteristic that these are very nearly the equations Newton 
would have written for an expanding gas acting under the influence of its own gravity. Note that 
in an expanding gas (in P = 1 models), one can redefine the Newtonian potential in such a way 
that equation (^) holds (see, e.g., Peebles, 1980; section II.7). The ubiquitous expression 

(50) 

can be interpreted in a simple Newtonian way. The first two terms are simply the acceleration 
(in the rest frame of the matter) due to the pressure gradient force, and the last term is the 
acceleration due to the divergence of the anisotropic part of the stress tensor. In the case of a 
universe filled with cold matter, the equations for V and A have exactly the same form as the 
Newtonian equations in expanding coordinates. This is why we refer to our treatment as the 
“quasi-Newtonian” formulation. 

This formulation is closest to working in the gauge known as the “velocity-orthogonal isotropic 
gauge” (Kodama & Sasaki, 1984) in which there is no residual gauge freedom. We will elaborate 
more on the physical meaning of these variables and the nature of this gauge in a future publication 
(Lyth & Schaefer, in preparation). In both the Aca and Va — V equations the term (^) appears. 
This term is in turn equal to 

which is effectively zero in a cold matter dominated universe. The quantity is the 

perturbation of the expansion rate due to the matter perturbations. 


1 


(1 -|- rc) a 


c, A -|- rcP-tell 

3 


k 


o(l -|- tc) 


c, A -|- tcT-tell 

3 
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Note that all of our evolution equations depend only on the perturbations in the distribution 
of particles. We do not need an explicit solution of the behavior of the metric perturbation as 
in the gauge dependent approach (Bond &: Szalay, 1983; Ma & Bertschinger, 1995). We also do 
not need an explicit solution of the potential equation (for $) as in Durrer & Straumann, (1988) 
or Stompor, (1994). All of the background cosmological information is carried in the equation 
for the expansion rate. To solve these equations then, the only remaining piece of information 
we need is the initial conditions. Before getting to these, we would like to derive a set of fluid 
equations for purely relativistic particles, as this would represent a great simplification over the 
equations in this section (for this component). 


3 Equations for Collisionless Relativistic Particles 


If the particles are effectively massless over the range of time of interest in the universe, so that 
u = g is at least a good approximation, we can integrate the equations for the cr^ to get equations 
for fluid variables. In this case Ta will be zero. 

We multiply the equations for at by v^, set q = v and then integrate over v to get the fluid 
variables. The first three fluid variables Aca, 14, and Ha have been defined previously. We repeat 
them again here for the case of relativistic particles for completeness. Using the subscript r for 
relativistic fluids: 

T^4 roo 

dvv^ao 


Aar — 


Ur = - 


rCrllr = 


H f 

Pr Jo 

3ri 


Ap. 

^4 po> 
- 
Pr d 0 

Tr = 0 


3 

/ dvv ui 
Jo 


dvv'^a2 


and we also have higher momentsf] which we call corresponding to ai for £ > 3: 


w. 


= 

* A 


rpA 
J. r- 


poo 

/ dvv^af, I > 3 
Jo 


(52) 

(53) 

(54) 

(55) 

(56) 


pr J 0 

If one were observing the radiation pattern from the origin of this coordinate system at time 
r, the would be the amplitudes of the coefficients of a multipole moment expansion of that 
pattern. 

The lowest moments equations are coupled by gravity to the other matter in the universe. 
The equations for these are 


4 4 d 

Ur T — PVr —- 

3 l + wa 


c, A + wT -uiB 


= 0 


(57) 


^Note that although the same symbol is used in Schaefer (1991), these are defined to be a factor of 3 larger 
than in Schaefer (1991). 
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(58) 



4 Initial Conditions: Solutions for Super—Horizon Perturbations 
In the Radiation Dominated Epoch 

To numerically integrate structure formation models, one requires an accurate set of solutions 
which may be used as initial conditions. To that end, we consider solutions from both of the 
“independent” types of perturbations, adiabatic and isocurvature, on super-horizon scales. Adi¬ 
abatic fluctuations typically arise from inflationary models and isocurvature modes often result 
from phase transitions in the early universe. In each type we will solve for the behavior of the 
growing mode of the perturbations well outside the horizon, deep within the radiation dominated 
era. These solutions will be the most relevant for starting numerical calculations. We will start 
first with the simpler adiabatic fluctuations. 

4.1 Adiabatic Perturbations 

In an adiabatic perturbation, the density fluctuations 6pa in in the different components (a, b, 
etc.) have amplitudes which are related by 

Spa _ Spb _ ... _ 

{l+Wa)pa {l+Wb)pb {l+w)p 

in order that the entropy per particle does not change. With our choice of gauge invariant density 
fluctuation variable, the same condition applies even when we are outside the horizon. Thus the 
adiabatic initial conditions are exactly as expected intuitively; 

Aca ^ Acfe ^ ^ A 

(l + Wa) (l + rt;;,) (l + w)’ 

for all components. We now proceed to identify the growing mode solutions since these are the 
only ones that will survive long after their formation. It will be assumed that we are beginning 
our simulation deep enough in the radiation dominated epoch that the temperature is much larger 


( 62 ) 


( 61 ) 
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than the mass of any hot components. Look first at eq. (|^). Because we are considering super¬ 
horizon scales, we can ignore all terms of order (A:r)A, an assumption which may be validated 
self-consistently. It is obvious that I 4 = IL is a consistent solution for all components: 


Va = Vb = --- = V. 


(63) 


Now consider eq. (^5|) , the equation for the growth of Aca- For each component, F^ and (c^ — Wa) 
are both negligible, and, since the temperature is much larger than the mass of any hot component, 
the factor l + tCa is a constant for each component as well. So, eq. (^) shows that the evolution of 
Aca/(1 + Wa) is the same for each component, preserving the adiabatic initial conditions. We are 
left with the greatly simplified problem of solving only one set of component equations. One can 
see from the equations for the evolution of the collisionless particles that we will need to retain 
the first three moments (A, V and 11). For the photons, due to their strong interaction with the 
baryons, and for the cold components, we only retain the first two moments. The total anisotropic 
pressure, H, which appears in the evolution equations for A and V can be found by solving for 
the anisotropic pressure of the collisionless components, Her- The total II is related to the the 
collisionless components by the ratio of the collisionless to total pressure, i.e. II = {per/p)^cr where 
Per is the pressure contribution from the collisionless components and p is the total pressure. It is 
left as a exercise for the reader to verify that the solutions for curvature perturbations have the 
following dependence on r: 


A{k,T) 

V{k,T) 

U{k,T) 


AHCi{k)k^T‘^ 

3 A{k, t) 

4 kr 

-^A{k,T) 


1 + 
1 + 


2pcr j ^ 

5p 

2pcrl“^ 


(64) 

(65) 

( 66 ) 


where a{k) is a stochastic variable which carries information about the probability and spectrum 
of the initial fluctuations, and A// is the amplitude of density fluctuation at Horizon crossing, 
i.e., when kr = 1. a{k) is usually assumed to be a Gaussian random variable [i.e. with random 
phase - see Bardeen, et al., 1986], and was inspired by the formalism of quantum field theoretic 
description of the quantum fluctuations during inflation (Guth &: Pi, 1982; Abbott k. Wise, 1984). 
When averaged over an ensemble of universes (or equivalently over many horizon volumes), a{k) 
has a mean of zero and a variance defined by 

{6t{k)a{k')) = {27r)^k^-^6^{k - k') (67) 


where the brackets denote the ensemble average and the * signifies the complex conjugate, n is 
the usual spectral index, and n = 1 defines the scale free (Harrison-Zeldovich) spectrum. The 
factors of 27r, which are not in Abbott & Wise, (1984), appear here because we have associated 
the l/(27r)^ factor with the wavenumber integral of the Fourier transform, instead of with the 
space integral transform. 

Note that Ah above differs from the more common specification of the amplitude at Hubble 
crossing when the wavenumber is equal to the Hubble length [ka/a = 1). Thus the amplitude 
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Ah is a factor of 4 larger than the parameter en of Abbott & Wise (1984). We have chosen to 
use the Horizon crossing specification so that Ah has the same definition in the radiation and 
matter dominated epochs. Using the more common Hubble crossing definition leads to the use of 
different amplitudes for radiation and matter dominated epochs (see e.g. Bardeen, Steinhardt, & 
Turner, 1983). In the strongly radiation dominated phase r is related to the scale factor by 

= 4.66 X lO^a, (68) 

where Hq is the present Hubble constant, and and are the present day matter and radiation 
(as if there were 3 types of relativistic neutrinos) energy densities. The above number assumes 
the COBE central value for the current photon temperature T = 2.726 K (Mather, 1994). 

The power spectrum P{k) of adiabatic perturbations is the Fourier transform of the two point 
density correlation function. Using our notation it can be simply computed from the following 
definition: 

{A*{k, T)A{k , r)) = {2'K)‘^P{k, T)6^{k — k ). (69) 

Our initial power spectrum at is then 

P{k,n) = {27rfAlk^Tt (70) 

The evolution equations in the previous sections can then be used to get the present day A{k, tq) 
and via eq (|6^) the present day power spectrum. 



4.2 Isocurvature Perturbations 


Isocurvature perturbations are distinctly different from the curvature type in that the variation 
is not in local energy density, but rather in local content. We specify the initial conditions as 
perturbations which violate eq. (^), in such a way as the total perturbation vanishes (A = 0). 
The relevant variable for parameterizing an isocurvature perturbation is then naturally 


Sab = 


Ar 


^cb 


l + Wa l + Wh 


(71) 


which measures the degree of non-adiabaticity. Note that in the adiabatic case there was only one 
set of initial conditions; here we see that isocurvature intial conditions allow for many different 
possibilities. Rather than discuss the system in general, it is more instructive to consider a specific 
type of isocurvature perturbation and see how to set up the isocurvature perturbations. It is 
typical to consider the type of perturbation where the perturbation in the matter is compensated 
by a perturbation in the radiation. For simplicity we will consider a case where the matter is 
baryons and the radiation is only photons. In this case the density perturbations will have the 
following form; 


1 + 3pf,/(4/9^) 

Pb Pbr 

Pr 1 T 3pft/(4pj.) 


(72) 

(73) 
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However, we are considering the initial conditions deep in the radiation dominated era, when we 
take {pb/Pr) 1- When solving the equations to find the growing mode, one must take care to 
keep only those terms which are of the same order in small quantities. For example, if we keep 
only the zeroeth order quantities in kr and {pb/Pr) in a universe with only baryons and radiation, 
the solution would be 

Aefe ~ Sbr (74) 

~ 0 A (75) 

V ^ O^Vr^Vb. (76) 


This is not accurate enough for numerical work, so we will find the leading order dependence to 
insure that roundoff error does not induce an unwanted adiabatic perturbation. In the specific case 
we are considering a small adiabatic perturbation is actually induced because the two components 
do not have the same sound speeds. In fact, it is generally true that adiabatic and isocurvature 
modes cannot be completely separated when different components have different sound speeds. 
(See, e.g., Kodama &: Sasaki, 1984). 

We now solve for the leading order terms in the initial conditions where we simply assume 
that A <C Aca. To make our universe more realistic, we will also include neutrinos which we 
assume to be distributed in the same manner as the photons {Sbr = Sbu)- 


A,r = A,, =- ^A,b (77) 

Pr + Pu 
Pb o 

pr Pv 

Because A is small, the term which drives the evolution of the fluctuations is the entropy pertur¬ 
bation, tcT, given in the relativistic epoch by 


T « 


Pb 


Pr T Pv 


Sbrt 


(78) 


The T terms in equations (^) and (^) then drive A and V. As was the case for adiabatic 
perturbations, only the first three moments are significant. One can easily verify that the solutions 
for isocurvature initial conditions are 


1 -1 


V{k,T) = ^kTT{k,T) 1 + ^^ 

n(fc, r) = —^^kTV{k,T) 

15 p 

A(/c, r) = — kTV{k,T) -n(fe,r) 


(79) 

(80) 
(81) 


Furthermore, it can be shown using eq.’s (^ &: ^) that 


Acr — Acu — F (82) 

v; = V^ = Vb = V. 
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As we will discuss in the next section, 14 and 14 are kept to be the same by the strong scattering 
of photons and baryons {via electrons). We can now clearly see from eq.’s (M) that our 
assumption that T induces total density perturbations A has been verified self consistently. 

We have avoided introducing a random variable in the isocurvature case because the fluc¬ 
tuations tend to come from non-gaussian sources. The probability distribution function is then 
specific to the particular source of the perturbations, and the accompanying formalism is perhaps 
best left unspecified. In the case of Gaussian intial fluctuations, the initial entropy perturbation 
Sab can be expressed using the stochastic variables used in the adiabatic case: Sab = Sa{k), where 
5* is a constant amplitude of the intial perturbation. The present day power spectrum can be 
directly computed using eq. (|6^) . 

5 Equations for Photons: Baryon Coupling and Temperature 
Anisotropies 

Observations of fluctuations in the temperature of the relic photons are an important diagnostic 
tool for studying theories of structure formation. The techniques for relating the measured cosmic 
microwave anisotropies to theoretical models often use the predicted multipole moment spectrum. 
We now want to specify how these moments are calculated using our formalism to facilitate 
comparisons of temperature anisotropies with estimates of the matter fluctuation amplitudes. 
Before doing this, however, we need to cover one more topic which is necessary for calculating the 
evolution of the photon distribution. 

Up to now we have only considered collisionless particles. When matter is ionized in the early 
universe, there is strong scattering between the photons and electrons. We discuss the treatment 
of this scattering in the next sub-section. 

5.1 Photons Coupled to Baryons 

The equations in section work well for photons when they are sufficiently decoupled from the 
baryons that they are collisionless. However in the early universe the temperature was high 
enough to keep normal (baryonic) matter fully ionized. The free electrons have a large cross 
section for interacting with the photons and the baryons, effectively coupling the baryons and 
photons tightly. As time progresses, the temperature cools off to the point where the electrons 
can combine with the protons (and nuclei) to form neutral atoms, which have a small cross section. 
We need equations to describe the coupling of the photons and baryons. The equations for this 
have been derived elsewhere (Kodama Sz Sasaki, 1984 - cf. Appendix E). Here we will just give the 
answer using the present formalism for the sake of completeness. For a more detailed explanation 
of some of the physical effects of on the cosmic photon distribution, we refer the reader to Hu &: 
Sugiyama (1995), who use the same density perturbation variables as in this work. 

The basic Liouville equation now has a collision term C{f), which depends on the distribution 
function /. The equation becomes 

m) = C{f). (83) 
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where C{f) is the scattering functional which must be worked out for the electromagnetic scat¬ 
tering of photons by electrons. To get the fluid equations, we must then integrate C{f) over 
momentum in the same way that we integrated the Boltzmann operator C. One can see that the 
equations are just the same as the collisionless equations, but now with an extra term added to 
describe the effect of scattering on the distribution. 

The equations for the photons are as follows. First, the equation for A^r is unchanged: 


T ——- 

3 l + wa 


ciA + wT -ten 

3 


= 0 


(84) 


The fluid velocity equation (as well as the higher moment equations), get scattering terms. 



o^r) + ^ f-') [A -I- 2wli] + -V 
3 2k \a/ a 


ane(7T{Vb - Vr) 


(85) 


where tt-e is the electron number density and or = 6.6524 x 10“^® cm^ is the Thompson scattering 
cross section. The product Uecrx is the photon collisional frequency. 

The higher moments go as 

Hr - 7(8Fr + 3n®) = -^aUefTT^r (86) 

5 ID 


and 




k 

2£ + 1 




{I +1)4^+^) 


— atlefTT’ 



(87) 


We have included the cos^ Q angular dependence in the differential scattering cross section. This 
is not explicitly evaluated in the formula in Appendix E of Kodama & Sasaki, but is included in 
their later work (Kodama &: Sasaki, 1986; Gouda &: Sasaki, 1986). 

The equations for the baryons can be obtained by knowing that the total energy momentum 
tensor is conserved. This means that the baryons behave like cold matter with a scattering term 
on the right hand side of the velocity equation. 


k 


4 Pj, 


[H - K] + -[14 - K] = -^_[c^A + wr- -?nn] + -^aUeaTiVr - Vb) 


a' ' 1 + rc ^ ^ ' 3 

while the density perturbation equation remains unchanged: 


3 Ph 


( 88 ) 


Acfo — —kVh + 


o ^ r 9 * ^ 

T—- [ci^ + wT 

l + wa 



(89) 


To see the effect of the coupling induced by the scattering terms it is instructive to look at the 
equations for the differences between the photon and baryon fluid variables. This is done most 
conveniently in terms of the variables Sbr and 14r. (Kodama &: Sasaki, 1984), where Sbr is the 
relative entropy perturbation as in eq. ( 0 ) 


St 


r 


^cb ^cr 

I + Wb I + Wr' 


(90) 
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and 


Vfyr — Vj} ^r- 

After a little algebra we arrive at the equations 


^br — kVfyi' 


and 


Vbr = - 


+ an^aT 1 + 


4 

3 pb 


Vbr - -{Vr - - 7 

a 4 




(91) 


(92) 


(93) 


The key to understanding the behavior of the photon baryon system is contained in the first term 
on the right hand side of eq. (93). If the collision frequency (neor) is much larger than the 
expansion rate {a/a?'), then the two fluids are tightly coupled. Note that in this regime small 
differences between the baryon and photon fluid velocities will be strongly damped out. As long 
as Vbr is very small, then the baryon density perturbation and the photon perturbation will be 
kept proportional [in adiabatic perturbation cases, Acfe = (3/4) Acr]. 

Through decoupling, it is numerically easier to integrate the equations for the photons and Sbr 
and Vbr- Well before decoupling we can obtain values for Vbr and 11^ by assuming “equilibrium” 
conditions, that is, by using the value of Vbr and 11^ for which Vbr = 0 and 11^ = 0 in equations 
( p^ & ^), respectively. We have found that the tight coupling limit (with these equilibrium 
solutions) is valid down to a temperature of about 6000 K below which one should integrate 
the full equations. After decoupling, we switch from integrating equations for Sbr and Vbr in 
favor of those for the baryon perturbations Acb and H- ^cb ( 14 ) can be recovered using a linear 
combination of Sbr ( 14 r) and /S.cr (W)- 


5.2 Evolution of Photon Perturbations in the Matter Dominated Epoch 

The equations introduced in the previous section are appropriate for determining the evolution of 
the photon perturbations through the epoch of recombination. However, to accurately track the 
photons until the present epoch requires the retention of an inordinate number of terms in the 
multipole expansion and is therefore numerically undesirable. Fortunately, following recombina¬ 
tion, it is possible to calculate an analytic solution by making a few reasonable approximations. 
Most significantly, the interactions between the photons and baryons can be ignored, and there¬ 
fore the collisionless Boltzmann equation may be solved. In the next sub-section, we will estimate 
the errors made by using this simplification. 

Since we are primarily interested in calculating the temperature fluctuation in the microwave 
background, it is convenient notationally to write our equations in the variable M = AT/T. 
Given AT/T = (l/4)Ap/p, Ai is related to the photon phase space distribution T by 

Ai = / Airdw^T, (94) 

4 Pr 4o 

where the multipole moments Aii are dehned similar to the a£ of T in eq. (^^, i.e. 

■t OO 

M = — Y.{21 + \)?P£{p)M£. (95) 

£=0 
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The variable M is more convenient to work with here, and it is related to the variables of section 
m as 


II 

0 

5 

7 

17 

21 

1 

4Acr 

(96) 

(47r)-^Mi = 

1 

-A-- 

(97) 

II 

(M 

— n,. 

12 

(98) 

II 

1 

-nw 

12 

(99) 


Cosmic microwave anisotropy observations can be expressed in terms of multipole moments de¬ 
rived from the spatially averaged temperature autocorrelation function (looking along two direc¬ 
tion vectors xi and ^ 2)5 

/ AT AT \ 

(^ — {xi) — {x2)j=^QlPl{xi-X2). (100) 


One also sees this series in terms of the Abbott &: Wise (1984) multipole moment coefficients 
which are related to the multipole moments by 


ai — . 


( 101 ) 


Note that the of used in Gorski, et ai, (1994) are a factor of l/(2£-|- 1) times a| as defined above. 
The theoretical multipole moments, Q'j, are found by integrating over all the Fourier components 
of Ail, i.e. 

1 fOO I 1 

( 102 ) 




However, while the Aif's are directly related to the observed quantities, it is more convenient for 
solving the Boltzmann equation to consider an alternate gauge invariant variable 0 (Kodama & 
Sasaki, 1986), defined by 

1 a. 


e = Ai- 

where and 4^ are the gauge-invariant potentials. T is given by 


--V 

k a 


+ $ = _3( y;n 


(103) 


(104) 


and <1> is defined in eq. (|4^ ). Using this definition of 0 along with the relation for the derivative 
of expressed in eq. (|5l|), the collisionless Boltzmann equation for the photons may be recast 
into 


0 -|- ik^Q = 


T - 4> 


(105) 


This equation can be solved by utilizing the integrating factor exp{ikfiT) producing the result 




T - 4> 


(106) 
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We designate the first term in this solution as the intrinsic term since it contains the intrinsic 
fluctuations as well as the Sachs-Wolfe effect. The second term is called the integrated Sachs- 
Wolfe (ISW) effect, and we will discuss it first. For the case of adiabatic perturbations in the 
matter dominated epoch, the potentials ‘h and T are effectively constant over the length scales 
of interest. One then can reasonably approximate the ISW effect by ignoring the oscillatory part 
of the integral because most of the change in the derivative occurs when A:(r — r') <C 1, i.e. the 
change occurs before the oscillations become significant. When this condition is violated, the SW 
effect is small relative to the intrinsic fluctuations and will not affect the approximation. The 
resulting ISW effect is then given by 

[($_^)_ (107) 


with At = T — Ti- The contribution to each multipole can be found by utilizing the Legendre 
expansion of the exponential 

CXD 

+ l)(-i)^j,(fcAr)P,(M), (108) 

0 

where ji is the spherical Bessel function. When considering isocurvature perturbations, the Sachs- 
Wolfe effect is small in proportion to the intrinsic fluctuations, so one can ignore the ISW effect 
altogether. 

Calculating the contribution from the intrinsic term is more difficult, as it involves finding the 
moments of products of more than 2 Legendre polynomials. These have been calculated by Bond 
& Efstathiou (1987) who showed that 


^ ^ 0*,(2r + l)CM.,™(-l)'-"^j;+zm2^(A:Ar), (109) 

l' m=0 


where the constant C is given by 




Q'—mCmQ—m 2/ T 2^ Am + 1 

cijf-i'—Yn + ‘IV — ‘Im + 1 


( 110 ) 


with 

_ (2n- 1)!! 

C-n — j 

n! 

The advantage to using this formulation should now be clear. One need only retain the number 
of terms significant after decoupling, which will be many fewer than are needed for an accurate 
solution today, and calculate a single summation rather than solve a large set of coupled differential 
equations. 



5.3 Adiabatic Sachs-Wolfe Multipole Moments 

A useful approximation of the multipole moments which can be evaluated analytically is the 
case of the Sachs-Wolfe effect for adiabatic perturbations. The simplest versions of inflation 
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predict Gaussian adiabatic fluctuations so this case is of particular interest theoretically (Abbott 
& Wise, 1984, Abbott & Schaefer, 1986; Fabbri, Lucchin, & Matarrese, 1987). For this reason 
this approximation has been used extensively in the analysis of COBE data (Wright, et al., 1992; 
Smoot, et al. 1992, Gorski, et al., 1994; Bennett et al., 1994). 

Our goal is to calculate the anisotropy in the photon distribution M at the present time tq. 
From eq. ( p.03| ), 

1 a. 


■M. 1 7-Q 01 ro 

During the matter dominated era a 


^ - V 
k a 


|ro¬ 


t^/tq and 


$ = 0 = 4', 


( 112 ) 


(113) 


which is true even for if the matter is made of free streaming particles like neutrinos as long as 
we stick to scales larger than the Jeans length k <C 0.02{mi^/l eV)^/^ h Mpc“^ (see e.g. Babu, et 
al., 1995), where rrii, is the neutrino mass. In this case the solution of the Boltzmann equation in 
eq. (|106|) takes a particularly simple form; 

©Iro = (114) 


where At = tq — t^, and is the time when the photons decouple from the baryons. 

On the largest scales, we can make the approximation M\t^ ~ 0. This is justified as follows. M 
is a linear combination of fluid variables Acr, Wj B^, and Hr . If /c <C then to leading order, 
Adlrd oc [kr^Y. On the other hand, the potential T oc {krd)^, so Al|rd ^ '!'• This approximation 
describes the redshifting of photons by the gravitational potentials of matter perturbations, which 
is often identified as the Sachs-Wolfe effect. Using this approximation. 
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During the matter dominated epoch. 


1 a. 


T = —6Ajtq;(A:), and —-U = —^/S.Ha{k). 
k a 

Using eq. (p^), then we find for i > 1: 

M(\ro = (-l)^+^87rA//Q!(fc)j>(A:Ar). 

For the multipole moment predictions we need the rms values of 

(MXk)M^(k')} = (27rflM^I^6(k-k') 

Using equation (p^) and plugging the resulting value of in eq. (|102|) we get 


^ F[2-n/2] r[£-n/2 + 5/2] 


(115) 
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In the special case of the scale free Harrison-Zeldovich spectrum (n = 1), this reduces to 


The above equation shows why the quantity i{l + 1)Q|/(2I' -|- 1), [called + l)C£/(47r) by some 
authors (e. 5 . Bond &: Efstathiou, 1987; Crittenden, et al., 1994; Hu &: Sugiyama, 1995; White, 
Scott, & Silk, 1994)] is often plotted versus i. If the Sachs Wolfe effect were the only source of 
anisotropy, this quantity would correspond to a horizontal line. Deviations from a straight line 
then illustrate the contributions from Afjrd and the integrated Sachs Wolfe effect. We will plot 
the results of a numerical calculation in section p. 


5.4 The Effects of Relic Electrons and Reionization 


We have so far ignored the effects of scattering on our results for the evolution of photon per¬ 
turbations in the decoupled epoch, but it may have significant consequences. Recombination is 
not complete and there is always present a small relic ionization. The observations of Gunn &: 
Peterson (1965) indicate reionization has occurred prior to a redshift of z ~ 5 (see Giallongo et 
al. 1995 for recent results). We would now like to consider these possibilities in our solutions and 
estimate their effects. 

The Boltzmann equation with interaction terms can be written in the following form: 


0 -|- ikfiQ + Rc—Q = 
a 


T - <I> 


T Rc . Rc Rc o T2(/^)) 

a 47 r a a an 


( 121 ) 


where Rc = a^neCTT/cL is the ratio of the interaction rate to the Hubble expansion rate and 0^ is 
related to 0 the same way as M.i is related to A4, i.e., through eq. (1^) . This leads to a modified 
integration factor which includes a real component, exp (^ikfin -|- , that tends to wipe 

out fluctuations on scales smaller than the scattering mean free path. Note that throughout 
this section we implicitly assume that oq = 1. Outside the horizon, as can be seen from eq.’s 

&: 1 ^ ), the photon perturbations grow like the dominant matter with scattering only affecting 
the growth of the H term; scattering effects become important on the sub-horizon scales. Well 
inside the horizon, the oscillating part of the integrating factor tends to eliminate the effects of 
the inhomogeneous terms in the Boltzmann equation, those appearing on the right-hand side of 
eq. (121). As a hrst approximation, it is reasonable merely to consider the homogeneous solution, 
that is to multiply the collisionless solution by the damping factor exp da^^ , where a* is 

defined as the greater of Oj and ahc, defined as the scale factor for which Hahc/k = 1. 

We now can use this approximate solution to estimate the effects of the relic ionization on 
our collisionless results. In the matter dominated epoch following recombination, we can write 
an expression for the damping factor Rc = .069(1 — Yp)Q.},hxea ~‘^^‘^where Yp = .23 is the Helium 
mass fraction and Xe is the Hydrogen ionization fraction. As a particular example, consider a 
universe with h = \l‘l and = .05 which leaves a relic ionization of approximately Xe ~ 6 x 10“^. 
Recombination occurs at a redshift of Zrec ~ 1250, so if one uses the collisionless solution beginning 
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at z = 300 the damping factor is .997, meaning that at most the error is less than .3%. One makes 
only a negligible error by neglecting the relic ionization. 

If the universe is reionized, we can also estimate the damping effect. Using the same model 
as the previous calculation, we assume that Xe = 1 for redshifts less than z ~ 5, consistent with 
Giallongo et al.(1995). We find a damping factor of .988 or a 1.2% reduction. 

6 Numerical Solution of the Equations 

We have implemented the equations presented here in two separate codes which agree very well 
with each other and with other calculations (Holtzman, 1989; Ma & Bertschinger, 1995). The first 
code is a fixed stepsize predictor corrector Haming type integrator which is extremely accurate 
when integrating oscillatory equations. The second is an adaptive stepsize Gear’s method implicit 
differentiation solver which offers a good balance between speed and accuracy. In the predictor 
corrector method, the stepsize is chosen for each specific k value, as the photon-baryon fluid 
oscillates with period 27r/k in conformal time. This integrator was reasonably fast. In the second 
method, we actually recast the equations here as differential equations where the variable was the 
scale factor a. During the matter dominated epoch a is not a linear function of r, so the stepsize 
needs to be changed continuously. This approach should of course only be taken using with a 
method with an adaptive stepsize. If we solve the perturbation equations as a function of a we 
need only specify the initial a, whereas with the fixed stepsize method we also have to calculate 
the initial value of r from eq. (^). Apart from the differences mentioned above the two codes 
are otherwise identical. 

The equations for evolution of the fluids and the collisionless particles are integrated explicitly 
to follow the evolution of density perturbations in the early universe. For the collisionless particles 
we integrate the at 15 momenta values and then use Gauss-Laguerre integration to get the 
associated fluid variables. We found the value of 15 momenta gave reasonably accurate results 
(error in A.cu was around 1 part in 10®). The universe we have modeled contains 1 flavor of 
massive neutrinos (collisionless particles), two flavors of massless neutrinos, photons, baryons, and, 
in the adiabatic case, cold dark matter. We treated the total fluid perturbations as independent 
variables from the separate components which was especially useful when dealing with isocurvature 
perturbations. The baryons and photons are treated using a tight coupling approximation until 
near recombination. The ionization fraction evolution is then approximated by interpolation to 
a table of separate numerical calculations (see Peebles, 1993). After z = 300, we no longer track 
the evolution of the photons and neutrinos as individual components. Keeping them requires 
the numerical evolution program to follow an exceedingly large number of oscillations of the 
relativistic fluids after that time when the wavelength is small (~ 1 Mpc). The net effect is to 
lose the influence of the anisotropic pressure in the equations governing the evolution of the total 
perturbations. However, since the contribution goes like tcH, it is not significant in the matter 
dominated epoch. 

We now present some results from the calculations for illustration. First we present the 
evolution of perturbations in an adiabatic cold plus hot dark matter model (G+HDM) model with 
5% baryons and 25 % hot dark matter (massive neutrinos) and a Hubble parameter h = 0.5. This 
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model was noted for its attractive large scale structure properties long before the COBE results 
were known (Schaefer, Shah &: Stecker, 1989; see also Holtzman, 1989). In figure 1 we show the 
evolution of the gauge-invariant fluid variables with time r for k = 0.01 /i/Mpc. In the three 
panels a), b) and c), we present the (component) perturbations in density, velocity, anisotropic 
pressure, and neutrino entropy. We also present the total density and velocity perturbations. The 
growth on this large length scale is not very exciting, as the massive neutrinos {rriu = 5.9 eV) 
are effectively cold on this scale. The adiabatic conditions are preserved until the perturbation 
crosses the horizon at r = 1/k. The photon-baryon fluid begins to oscillate, but soon experience 
decoupling at r ~ 220 Mpc. The photons then, like the massless neutrinos, simply free stream 
after that time. The entropy in the neutrino distribution grows but never becomes significant, 
especially when we consider that it only enters the through the combination where Wi, 

decreases rapidly after the neutrinos become non-relativistic. We stop following the evolution of 
the relativistic components (photons and massless neutrinos) at r ~ 700 Mpc, which corresponds 
to a redshift of z 300. 

In figure 2 we present similar results for k = 1 h/Mpc. Here we start with the same initial 
amplitude as in figure 1: A(rj = 0.06 Mpc) = 0.01. The evolution is more interesting. First of 
all we see in panel a) that the CDM component experiences growth which is damped (relative to 
figure 1) by the effects of photons and neutrinos during the radiation dominated epoch and by 
the free streaming of the massive neutrinos in the matter dominated epoch. The time of radiation 
matter equality is Teq ~ 60 Mpc. The photon-baryon system undergoes acoustic oscillations 
maintaining the adiabatic condition A^b = (3/4) Ac,, until decoupling. Note that the velocity and 
density perturbations are 90° out of phase, as the particles flow in and out of the perturbation, 
the density contrast must lag the velocity. As the photons and baryons begin to decouple, they 
undergo viscous (Silk) damping. After decoupling, the baryons fall into the potential wells of the 
CDM component. The massive neutrinos free- stream, but slowly settle into the CDM potential 
wells as the neutrino momenta redshift. 

We next consider a case with isocurvature perturbations, the isocurvature hot dark matter 
(IHDM) model. We begin the model with the initial conditions as specified in eqs. (77 - 32). 
Again we plot the results in figure 3 for k = 0.01 h Mpc“^ and in figure 4 for k = 1 h Mpc“^. Here 
there is no CDM component and with h = 0.5 we have Di, = 1 with a neutrino of mass rriy = 23 
eV. In panel c) we plot the total entropy T, because of its role in driving the density perturbation 
growth. In figures 3 and 4 we see that A^b is approximately constant until the induced adiabatic 
curvature perturbation A becomes of the same order. Outside the horizon the neutrino and 
photon perturbations scale with the ratio of the matter density over the radiation density, while 
the total perturbation is down by a factor as in eq. ^ The velocities and anisotropic 
pressure perturbations scale with the amplitude of the induced adiabatic perturbations. 

In figure 4, again we see that A^b remains constant until the induced adiabatic perturbation 
amplitude becomes comparable. Now we also see that despite the different amplitudes of Acb and 
Acr the velocity perturbations oscillate together in the same way as in an adiabatic perturbation. 
We also see the viscous damping effect, and as in figure 2 the massive neutrinos gradually fall 
into the potential wells set up by the baryons. 

As an example of using our formalism for calculating temperature anisotropies, we also calcu- 
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late the multipole moments for Gaussian n = 1 spectra in an adiabatic CDM (figure 5) and 
an isocurvature HDM universe (figure 6). The multipole moments in figure 5 show nearly a hori¬ 
zontal line for the low as expected from the Sachs-Wolfe approximation of section 5.2. We have 
normalized the moments to COBE, using the quadrupole {Q2) with the value Qrms-PS = 20 ^ 

K. 

7 Conclusions 

We have presented a simple formalism for following the evolution of density perturbations in 
the early universe. Our treatment includes the important case of collisionless relic particles like 
massive neutrinos, which cannot be accurately described as a fluid. The density perturbations 
variables have properties which make their evolution appear to be Newtonian with an expanding 
space. This makes it much easier to understand the physics of the evolution, and we find this 
formalism to be very attractive for this feature alone. We have demonstrated the use of these 
equations for transfer functions and temperature anisotropies in both adiabatic and isocurvature 
models. 

We have shown examples here which demonstrate how the formalism works. The equations 
here can be easily modified to incorporate other cases, as we have shown elsewhere (de Laix, 
Scherrer, & Schaefer, 1995). The initial conditions are particularly easy to specify here (especially 
in isocurvature models). We have presented this new formalism as a working man’s set of equations 
which do not require a deep understanding of general relativity to comprehend the behavior the 
perturbations under the action of gravity. 

Lastly, if one is interested in approximate transfer functions we have presented a modified set 
of equations (in the Appendix) which can be integrated extremely quickly on a computer. 

Acknowledgements A. d. L. wishes to thank the DOE for support under the grant DE- 
AC02-76ER01545. R. S. wishes to acknowledge support for this research under NASA grant 
NAG5-2646. 

8 References 

Abbott, L.F., & Schaefer, R.K., 1986, ApJ, 308, 546 
Abbott, L. F., & Wise, M., 1984, ApJ, 282, L47 

Babu, K. S., Schaefer, R.K., & Shah, Q., 1995, Bartol Preprint BA-95-19 
Bardeen, J. M., 1980, Phys. Rev. D, 22, 1882 

Bardeen, J. M., Bond, J.R., Kaiser, N., & Szalay, A. S., 1986, ApJ, 304 , 15 
Bardeen, J. M., Steinhardt, P., & Turner, M.S., 1983, Phys. Rev. D, 28, 679 
Bennett, C. L., et al, 1994, ApJ, 436, 423 


25 



Bond, J. R., &: Efstathiou, G., 1987, MNRAS, 226, 655 

Bond, J. R., & Szalay, A.S., 1983, ApJ, 274, 443 

Crittenden, R. et al, 1994, Phys. Rev. Lett., 71, 324 

de Laix, A. A., Scherrer, R.J., &: Schaefer, R. K., 1995, ApJ, (in press) 

Durrer, R., 1989, A&A, 208, 1 

Durrer, R., &: Straumann, N., 1988, Helv. Phys. Acta, 61, 1027 

Ehlers, J., 1971, Proc. XLVII International School Enrico Eermi, ed B.K. Sachs, (Amster¬ 
dam; North Holland) 

Fabbri, R., Lucchin, F., & Matarrese, S., 1987, ApJ, 315, 1 

Giallongo, E., D’Odorico, S., Fontana, A., Savaglio, S., Cristiani, S., & Molaro, P., 1995, 
astro-ph /950300^ 

Gorski, K., et al, 1995, ApJ, 370, L5 

Gouda, N. Sz Sasaki, M., 1986, Prog. Theor. Phys., 76, 1016 

Gunn, J. E., & Peterson, B. A., 1965, Astrop. J., 142, 1663 

Guth, A. &: Pi, S.-Y., 1982, Phys. Lett., 116B, 335 

Holtzman, J. A., 1989, ApJS, 71, 1 

Hu, W. & Sugiyama, N., 1995, Phys. Rev. D, 51, 2599 

Kodama, H., &: Sasaki, M., 1984, Prog. Theor. Phys., SuppL, 78, 1 

Kodama, H., &: Sasaki, M., 1986, Int. J. Mod. Phys., A 1, 265 

Lyth, D.H., & Liddle, A. R., 1995, Astrop. Lett. & Comm, (in press) 

Lyth, D. H. &: Bruni, M., 1994, Phys. Lett., B323, 118 

Lindquist, R. W., 1966, Ann. Phys., 37, 487 

Ma, C.-P. & Bertschinger, E., 1995, ApJ, (in press) 

Mather, J.C. et al, 1994, ApJ, 420, 439 

Mukhanov, V.F., Feldman, H.A., & Brandenberger, R.H. 1992, Phys. Rep., 215, 203 
Peebles, P. J. E., 1973, ApJ, 180, 1 


26 





Peebles, P. J. E., 1980, The Large Scale Structure of The Universe. (Princeton University 
Press: New Jersey) 

Peebles, P. J. E., 1993, Principles of Physical Cosmology. (Princeton University Press: New 
Jersey) 

Pogosyan, D. Yu., & Starobinsky, A. A., 1994, MNRAS, 265, 507 

Rebhan, A. K., & Schwartz, D. J., 1994, Phys. Rev. D, 50, 2541 

Schaefer, R. K., 1991, Int. J. Mod. Phys., 6, 2075 

Schaefer, R. K., & Shah, Q., 1994, Phys. Rev. D, 49, 4990 

Schaefer, R. K., Shah, Q., & Stecker, P.W., 1989, ApJ, 347, 575 

Smoot, G.E., et al, 1992, ApJ, 396, LI 

Stewart, J. M., 1972, ApJ, 176, 323 

Stompor, R., 1994, A&: A, 287, 693 

Valdarnini, R., & Bonometto, S. A., 1985, A&A, 146, 235 

White, M., Scott, D. & Silk, J, 1994, AARA, 32, 319 

Wright, E. L., et al., 1992, ApJ, 396, L13 


9 Figure Captions 

1. A plot of temporal evolution (in conformal time r) of perturbations in a) density b) velocity 
and c) anisotropic pressure for the various components and neutrino entropy in an adiabatic 
cold + hot dark matter model (with = 0.25). The perturbation wavenumber is A: = 0.01 h 
Mpc, and we are using h = 0.5. We stop following the relativistic neutrino and photon 
perturbations after a redshift of ~ 300. 

2. Similar to figure 1, using the same model, but for k = I h Mpc. Here we can see the effects 
of viscous damping of the photon-baryon fluid, and the gradual falUin of the neutrinos (and 
the baryons after decoupling) into the CDM potential wells. 

3. Same as figure 1, but here we consider an isocurvature HDM model, which has no CDM 
component. The baryon density perturbation is compensated for by the perturbation in the 
relativistic (photons and neutrinos) components. The rapid rise of the induced adiabatic 
perturbation A can be seen. 

4. Same as figure 3, but for k = 1 h Mpc. As in figure 2 we can see the various subhorizon 
evolutionary processes taking place. 


27 



5. Multipole moments of the temperature anisotropy as a function of in an re = 1, h = 0.5, 
CDM universe, normalized to Q 2 = 20 fi K/Tq. If the anisotropy were only due to the Sachs- 
Wolfe effect, the multipole moments would fit a horizontal line up to the £ corresponding 
the to horizon size at matter domination. 

6. Same as figure 5, but for the isocurvature HDM universe. Here initial perturbations in the 
photon distribution dominate the anisotropy. The moments are normalized to Q 2 = 20 ^ 
K/To. 

7. Comparison of transfer functions calculated using the full treatment solid lines and the 

approximation described in the appendix dotted lines for = 1 — VLcdm — = 0, 0.1, 

0.2, 0.3, 0.4 and 0.5, using h = 0.5 and adiabatic initial conditions. 


A A Quick Method for Getting Transfer Functions 


Setting up and running the codes to integrate the equations presented here can be time consuming 
in terms of man-hours and CPU time. If one only needs an accuracy of about 10% in the mass 
fluctuation, the procedure described above can be replaced with a very simple set of equations. 
These equations have been discussed previously (Schaefer, 1991). Here we will summarize that 
method and compare it to the exact calculations of the current work. 

The method, based on the “Grad” approximation of kinetic theory relies on 2 simplifications. 
First, for the collisionless particle components of the universe, we ignore all the higher moments 
where {£ > 3). This means we only consider Ac,., U-, and H^, and their equations, and set 
H^ ' = 0, for £ > 3. This simplification, when used for the massless relics, is well known and has 
been exploited before {e.g., Kodama & Sasaki, 1986). However, in order to use it effectively for 
the massive neutrinos, one needs a second approximation. 

The reason for this is that the equation for the massive neutrino lii, is problematic. As de¬ 
scribed in Schaefer, (1991), one can derive the fluid equations from the sigma equations. However, 
the equations involve two different momentum space integrals of each ai (see eq. The cor¬ 
responding fluid variables for these two momentum space integrals for no are the density and 
pressure perturbations. The pressure perturbation shows up in the T term in the equations pre¬ 
sented here. If the particles are relativistic so that the energy and momentum are equal, these two 
integrals are identical, and the problem disappears. Schaefer (1991) advocated approximating the 
two integrals as 



av — a£, 
Q 


{£ = 0 , 2 ). 


( 122 ) 


Using this approximation, we find T = 0 and we can derive an equation for H: 


(=f)' 

The other two equations are 

Ac^ Y 

1 + rfjy/ 


-3-{cl - Wu) + 3fc(l -F Wu)‘^Vu. 

a \ ct J 5 

= -(c^A-FrcT-^wH), 

1 + w a 3 


(123) 


(124) 
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(c^A + uiF 

~{^iy^cu ~ 



( 125 ) 


These equations are certainly not adequate for calculating the temperature anisotropy, because 
we have made a gross error in treating the photons. However, the density transfer functions 
are reasonably accurate. In Figure © we show a comparison of transfer functions using the 
approximations in this appendix and the full equations. Over the interesting range fc = 0 — 
1.0 h/Mpc, the error is less than 10%. 
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Figure 1. k= 0.01 h Mpc — 
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Figure 2 . k= 1 h Mpc — Qi, 
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Figure 4. k= 1 h Mpc 
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Figure 5. Multipole moments (CDM) 
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